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We propose two models for the interpolation between RGB images based on 
the dynamic optimal transport model of Benamou and Brenier [8] . While the 
application of dynamic optimal transport and its extensions to unbalanced 
transform were examined for gray-values images in various papers, this is the 
first attempt to generalize the idea to color images. The nontrivial task to 
incorporate color into the model is tackled by considering RGB images as 
three-dimensional arrays, where the transport in the RGB direction is per¬ 
formed in a periodic way. Following the approach of Papadakis et al. [35] for 
gray-value images we propose two discrete variational models, a constrained 
and a penalized one which can also handle unbalanced transport. We show 
that a minimizer of our discrete model exists, but it is not unique for some 
special initial/final images. For minimizing the resulting functionals we apply 
a primal-dual algorithm. One step of this algorithm requires the solution of 
a four-dimensional discretized Poisson equation with various boundary con¬ 
ditions in each dimension. For instance, for the penalized approach we have 
simultaneously zero, mirror and periodic boundary conditions. The solution 
can be computed efficiently using fast Sin-I, Cos-II and Fourier transforms. 
Numerical examples demonstrate the meaningfulness of our model. 


1. Introduction 

Color image processing is much more challenging compared to gray-value image process¬ 
ing and usually, approaches for gray-value images do not generalize in a straightforward 
way to color images. For example, one-dimensional histograms are a very simple, but 
powerful tool in gray-value image processing, while it is in general difficult to exploit 
histograms of color images. In particular, there exist several possibilities to represent 
color images m- In this paper we consider the interpolation between two color images 
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Figure 1: Left: RGB image of size 4x2. Right: Image as three-dimensional 4x2x3 
array. The “mass” is the sum of the RGB values of all pixels. 


in the RGB space (see Figure]^ rnotivated by the fluid mechanics formulation of 

dynamic optimal transport by Benamou and Brenier |8] and the recent approaches of 
Papadakis et al. [35] and [33] for gray-value images. In these works gray-value images 
are interpreted as two-dimensional, finitely supported density functions /o and fi of 
absolutely continuous probability measures /io and /ii (i.e. = J^/^dx, i = 0,1). 

In particular, we have Jj^ 2 /odx = /i dx = 1. Therewith, intermediate images are 
obtained as the densities ft of the geodesic path d/i^ = ft dx with respect to the Wasser- 
stein distance between /io and fii. 

In this paper, we extend the transport model to discrete RGB color images. The incorpo¬ 
ration of color into the above approach appears to be a non trivial task and this paper is 
a first proposal in this direction. We consider Ni x N 2 RGB images as three-dimensional 
arrays in M^iW2,3^ where the third direction is the “RGB direction” that contains the 
color information; for an illustration see Figure Particular attention has to be paid 
to this direction and its boundary conditions. We propose to use periodic boundary 
conditions, which is motivated as follows: assume we are given two color pixels /o and 
fi G Using mirror (Neumann) boundary conditions in the third dimension, the 

transport of, e.g., a red pixel /o = (1,0,0) into a blue one fi = (0, 0,1) goes over (0,1, 0) 
(green), see Figure]^ (middle), which is not intuitive from the viewpoint of human color 
perception. Furthermore, it implies that the transport depends on the order of the three 
color channels, which is clearly not desirable. As a remedy, we suggest to use periodic 
boundary conditions in the color dimension. In case of a red and a blue pixel, it yields a 
transport via violet, which is also what one would expect, compare Figure]^ (right) and 
see also Figures and [^ 

In the following, we propose two variational models for the transport of color images. 
To handle also the case of unequal masses the mass conservation constraint is relaxed. 
Our first model contains the continuity equation as a constraint. It turns out that it 
generalizes the technique which was proposed in [35] for gray-value images. The second 
model penalizes the continuity equation, similar as it was also considered for (continu¬ 
ous) gray-value images in [33]. Other approaches to unbalanced optimal transport can 
e.g. be found in laiiHiEgEH]. 

The interpolation based on an optimal transport model for a special class of images, 
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Figure 2: Left: RGB color cube. Middle/Right: Color transfer between (1,0,0) (red) 
and (0, 0,1) (blue) for mirror (middle) and periodic (right) boundary conditions 
visualized in the RGB cube. 


namely so called microtextures, has been addressed by Rabin et al. in [37l |4T]. The 
authors show that microtextures can be well modeled as a realization of a Gaussian 
random field. In this case, theoretical results guarantee that the intermediate measures 
jit are Gaussian as well and they can be stated explicitly in terms of the means and 
covariance matrices learnt from the given images /o and fi. The idea can be generalized 
for interpolating between more than two microtextures by using barycentric coordinates. 
The approach fails for some special microtexture settings which usually do not appear 
in practice. 

Color interpolation between images of the same shape can be realized by switching to the 
HSV or HSI space. Then only the hue component has to be transferred using, e.g. by a 
dynamic extension of the model for the transfer of cyclic measures (periodic histograms) 
of Delon et al. [251 SO] . For an example we refer to m- This idea is closely related to the 
affine model for color image enhancement in [M]. However, these approaches transfer 
only the color and leave the original edge structure of the image untouched. 

Finally we mention that the interpolation between images can also be tackled by other 
sophisticated techniques such as metamorphoses, see |19]. These approaches are beyond 
the scope of this paper. A combination of the optimal dynamic transport model with 
the metamorphosis approach was proposed in [33] . 

The outline of our paper is as follows: in Section we recall basic results from the 
theory of optimal transport. At this point, we deal with general p G (1,2] instead of 
just p = 2 as in [35]. We propose discrete dynamic transport models in Section]^ 
Here we prefer to give a matrix-vector notation of the problem to make it more intu¬ 
itive from an linear algebra point of view. We prove the existence of a minimizer and 
show that there are special settings where the minimizer is not unique. In Section 
we solve the resulting minimization problems by primal-dual minimization algorithms. 
It turns out that one step of the algorithm requires the solution of a four-dimensional 
Poisson equation which includes various boundary conditions and can be handled by 
fast trigonometric transforms. Another step involves to find the positive root of a poly- 
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nomial of degree 2g — 1, where ^ ^ = 1 ^ind p G (1,2]. For this task we propose 

to use Newton’s algorithm and determine an appropriate starting point to ensure its 
quadratic convergence. Section shows numerical results, some of which were also re¬ 
ported at the SampTA conference 2015 m More examples can be found on our website 
http://www.mathematik.uni-kl.de/imagepro/members/laus/color-OT, Finally, Section 
contains conclusions and ideas for future work. In particular, additional priors may be 
used to improve the dynamic transport, e.g., a total variation prior to avoid smearing 
effects. The Appendix reviews the diagonalization of certain discrete Laplace operators, 
and provides basic rules for tensor product computations. Further it contains some 
technical proofs. 


2. Dynamic Optimal Transport 


In this section we brieffy review some basic facts on the theory of optimal transport. For 
further details we refer to, e.g. nmsuso]. 

Let be the space of probability measures on and Tp(M^), p G [l,oc) the 

Wasserstein space of measures having finite p-th moments 

Vp{R^) := G V{R^) : [ \x\Pdp{x) < +oc 

I JRd 

For G let n(/io,/ii) be the set of all probability measures on x 

whose marginals are equal to po and pi. Therewith, the optimal transport problem 
(Kantorovich problem) reads as 

argmin / \x — y\^ d'K{x^y). 

7reIl{fio,fii) Jr^ 


One can show that for p G [l,oc) a minimizer exists, which is uniquely determined for 
p > 1 and also called optimal transport plan. In the special case of the one-dimensional 
optimal transport problem, if the measure po is non-atomic, the optimal transport plan 
is the same for all p G (l,oc) and can be stated explicitly in terms of the cumulative 
density functions of the involved measures. The minimal value 

Wp{fio,ni) := { min / |a: - d7r(a:, y) 

\7TeU{flo,fll) J^d 

defines a distance on 'Pp(M^), the so-called Wasserstein distanee. 

Wasserstein spaces (7^p(M^), IFp(M^)) are geodesic spaces. In particular, there exists for 
any po^ pi G Vp{W^) a geodesic 7 : [0,1] ^ Vp{W^) with 7 ( 0 ) = po and 7 ( 1 ) = pi. For 
interpolating our images we ask for pt — 7 (^), t G [0,1]. 

At least theoretically there are several ways to compute pt. If the optimal transport plan 
TT is known, then pt = yields the geodesic path, where xM^ ^ M^, 

Ct{x^ y) = (1 —t)x + ty is the linear interpolation map, see further [43]. This requires the 
knowledge of the optimal transport plan tt and of . At the moment there are efficient 
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ways for computing the optimal transport plan tt only in special cases, in particular in 
the one-dimensional case by an ordering procedure and for Gaussian distributions in the 
case p — 2 using expectation and covariance matrix. For p — 2 one can also use the 
fact that TT is indeed induced by a transport map T: ^ i.e., tt = (id,T)^/io, 

which can be written as T = V'0, where ^ fulfills the Monge-Ampere equation, see jl5| . 
However, this second order nonlinear elliptic PDF is numerically hard to solve and so 
far, only some special cases were considered HKniEniES]. Other numerical techniques 
to compute optimal transport plans have been proposed, e.g., in [ammaiiii. Another 
approach consists in relaxing the condition of minimizing a Wasserstein distance by using 
instead an entropy regularized Wasserstein distance. Such distances can be computed 
more efficiently by the Sinkhorn algorithm and were applied within a barycentric ap¬ 
proach by Cuturi et al. [211122]. 

The approach in this paper was inspired by the one of Benamou and Brenier in [8]. 
It involves the velocity field v: [0,1] x ^ of the geodesic curve joining po and 
pi. This velocity field v(t, •) has constant speed ||v(t, = Wp{po^ pi). It can be 

shown that it minimizes the energy functional 

Jo jRd. p 

and fulfills the continuity equation 


dtpt + Va, • •)) = 0, 

where we say that t is a measure-valued solution of the continuity equation if for 

all compactly supported 0 G G^((0,1) x and T G (0,1) the relation 

[ [ dt(l){x,t) + (v(x,t), V^(/)(x,t))d/it(x)dt = 0 

Jo JRd 

holds true. For more details we refer to [13] . 

Assuming the measures po and pi to be absolutely continuous with respect to the 
Lebesgue measure, i.e. dpi = /zdx, z = 0,1, theoretical results (see for instance m 
Theorem 8.7]) guarantee that the same holds true for dpt = /tdx, where ft can be 
obtained as the minimizer over v, / of 

^p(v,/)= [ [ -\v{x,t)\Pf{x,t) dxdt (1) 

JO JR'i P 

subject to the continuity equation 

dtf{x,t) + Vx ■ {v{x,t)f{x,t)) = 0, 

/(0,-) = /o, /(I, •) = /!, 

where we suppose supp/(t, •) C [0,1]*^ with appropriate (spatial) boundary con- 

ditions. Unfortunately, the energy functional 0 is not convex in / and v. As a 
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remedy, Benamou and Brenier suggested in the case p — 2 d. change of variables 
(/^v) (/?/v) = This idea can be generalized to p ^ (l,oo), see [13], which 

results in the functional 


fl 

Jo Jr^ 


Jp(rn{x^t)^ /(x, t)) dx dt, 


( 2 ) 


where x M ^ M U {+oc} is defined as 


y) 


if7/>0 

pyP-l II V > '-'5 

0 if (x,y) = (0,0), 

+OC otherwise 


(3) 


and I • I denotes the Euclidean norm. This functional has to be minimized subject to the 
continuity equation 


dtf{x,t) + \/x-m{x,t) = 0, (4) 

/(O,-) =/o, /(I,-) =/i, (5) 

equipped with appropriate spatial boundary conditions. 

Remark 1. The function J^: x M ^ M U {+oc} defined in © is the perspective 

function offj{s) = i.e., Jp{x^y) — yf; For properties of perspective functions 

see, e.g., [23]. In particular, since is convex for p G (l,oo), its perspective Jp{x,y) is 
also convex. Further, Jp{x,y) is lower semi-continuous and positively homogeneous, i.e. 
Jp{\x, \y) = \ Jp{x, y) for all A > 0. 

3 . Discrete Transport Models 

In practice we are dealing with discrete images whose pixel values are given on a rectan¬ 
gular grid. To get a discrete version of the minimization problem we have to discretize 
both the integration operator in Q by a quadrature rule and the differential operators 
in the continuity equation The discretization of the continuity equation requires 
the evaluation of discrete “partial derivatives” in time as well as in space. In order to 
avoid solutions suffering from the well known checkerboard-effect (see for instance [36] ) 
we adopt the idea of a staggered grid as in m , see Figure 

The differential operators in space and time are discretized by forward differences, and 
depending on the boundary conditions this results in the use of difference matrices of 
the form 
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Figure 3: Staggered grid for the discretization of the dynamic optimal transport problem, 
where N — P — A. For periodic boundary conditions, the boundary values for 
m are equal, while they are zero for Neumann boundary conditions. 



For the integration we apply a simple midpoint rule. To handle this part, we introduce 
the averaging/interpolation matrices 
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Discretization for one spatial dimension + time: In the following, we derive the 
discretization of Q for one spatial direction, i.e., for the transport of signals. The prob¬ 
lem can be formulated in a simple matrix-vector form using tensor products of matrices 
which makes it rather intuitive from the linear algebra point of view. Moreover, it will 
be helpful for deriving the fast trigonometric transforms which will play a role within our 
algorithm. The generalization of our approach to higher dimensions is straightforward 
and can be found in Appendix 0 

We want to organize the transport between two given one-dimensional, nonnegative 
discrete signals 

/o := fi := • 

We are looking for the intermediate signals ft for t — p, k — 1 ,..., P — 1. Using the 
notation ft{x) = /(x, t), we want to find ^ For m we have to 

take the boundary conditions into account. In the case of mirror boundary conditions, 
there is no flow over the boundary and since m = /v the value of m is zero at the 
boundary at each time. In the periodic case both boundaries of m coincide. The values 

k— - 

of m are taken at the cell faces ..., — 1 and time —p^, A: = 1,..., P, i.e., we 

/ . 

are looking for ( m(^, —p^) ) G where 

\ / j=K^k=l 

( 1 mirror boundary, 

[ 0 periodic boundary. 

The midpoints for the quadrature rule are computed by averaging the neighboring two 
values of m and /, respectively. To give a sound matrix-vector notation of the discrete 
minimization problem we reorder m and / columnwise into vectors vec(/) G 
and vec(m) G which we again denote by / and m. For the vec operator 

in connection with the tensor product 0 of matrices we refer to Appendix More 
specifically, let In G be the identity matrix and set 

Pf 


^ I Ip ® Sft mirror boundary. 

Ip ® periodic boundary. 


Pm := 


Ip ® mirror boundary. 

Ip ® (P^^)^ periodic boundary. 


Finally, we introduce the vectors 


1 


o(/o , / ■= P{-fS^O,p) , 






where we denote by 0 (and 1) arrays of appropriate size with entries 0 (and 1). They 
are used to guarantee that the boundary conditions are fulfilled. Now the continuity 
equation Q together with the boundary conditions Q for / can be reformulated as 
requirement that (m, /) has to lie within the hyperplane 


Co 



r}- 


( 6 ) 


We will see in Proposition that ^4^4^ is rank one deficient. Since further l^A = 0 , 
we conclude that the under-determined linear system in Q has a solution if and only if 
1^ f~ = 0, i.e., if and only if /o and fi have the same mass 




(7) 


This resembles the fact that dynamic optimal transport is performed between prob¬ 
ability measures. The interpretation of a color image as a probability density function 
has a major drawback; to represent a valid density, the sum of all RGB pixel values of 
the given images, i.e., the sum of the image intensity values, has to be one (or at least 
equal). Therefore, we consider more general the set 


C := argmin||(L>m|-Df) (7^) - / 111- (8) 

(mj) \J / 

Note that the boundary conditions (§ for / are preserved, while the mass conservation 
0 is no longer required. Clearly, if 0 holds true, then C coincides with Cq. Let lq 
denote the indicator function of C defined by 


I'Cix) 


0 if X G C, 
-hoo if X 0 C. 


For p G (1, 2], we consider the following transport problem: 


Constrained Transport Problem: 

argmin£'(m,/) := IVp(S'mm, S'f/+/+)||i + /)• (9) 

(mJ) 

Here, the application of Jp is meant componentwise and the summation over its (non¬ 
negative) components is addressed by the £i-norm. The interpolation operators S^n and 
Si arise from the midpoint rule for computing the integral. 

We can further relax the relaxed continuity assumption (m, /) G C by replacing it by 


where r > tq := niin(^ ^ — /~|| 2 - For t = tq we have again problem 

0 . Since there is a correspondence between the solutions of such constrained problems 
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with parameter r and the penalized problem with a corresponding parameter A, see 
13 lia SHI, we prefer to consider the following penalized problem with regularization 
parameter A > 0: 


Penalized Transport Problem: 


argmin Ex{m^ f) 
{mj) 


+ f+)h + X\\{D^\Di) (j^-r\\l ( 10 ) 


Note that also for our penalized model the boundary conditions for / still hold true. 
In general, both models Q and (10) do not guarantee that the values of / stay within 
the RGB cube during the transport. This is not a specific problem for color images but 
can appear for gray-value images as well (gamut problem). A usual way out is a final 
backprojection onto the image range. An alternative in the constrained model is a simple 


modification of the constraining set in ^ towards C := argmin^j^jQ \\{D^\Df) 



/“II2. This leads to inner iterations of the Poisson solver and a projection onto the cube 
in the subsequent Algorithm 1. In the penalized problem, the term ^[o,i]3 could be added. 
However, we observed in all our numerical experiments only very small violations of the 
range constraint, which are most likely caused by numerical reasons. 

A penalized model for the continuous setting and gray-value images was examined in 
[33] . For recent papers on unbalanced transport we refer to [ZKISlEniEHl. 

To show the existence of a solution of the discrete transport problems we use the concept 
of asymptotically level stable functions. As usual, for a function F: ^ M U {+00} 

and fi > iuG T(x), the level sets are defined by 


lev(F, /i) := {x G : F{x) < /i}. 

By Foo we denote the asymptotic (or recession) function of F which according to [23], 
see also [U Theorem 2.5.1], can be computed by 

\ T. . r- Fitx') 

Foo{x) = hmmf-. 

x' t 

t^oo 

The following definition of asymptotically level stable functions is taken from [U p. 94]: 
a proper and lower semicontinuous function F: ^ MU {-hoo} is said to be asymptot¬ 
ically level stable if for each p > 0, each real-valued, bounded sequence and each 

sequence {xk}k satisfying 

Xk G lev(F,/i^), ^ +OC, ^- > X ^ ker(Foo), (H) 

W^kh 

there exists ko such that 

Xk — px ^ lev(F, pk) foi* all k > ko. 
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If for each real-valued, bounded sequence {iik}k there exists no sequence {xk}k satis¬ 
fying then F is automatically asymptotically level stable. In particular, coercive 
functions are asymptotically level stable. It was originally exhibited in [5] (without 
the notion of asymptotically level stable functions) that any asymptotically level stable 
function F with inf F > —oo has a global minimizer. A proof was also given in [H 
Corollary 3.4.2]. With these preliminaries we can prove the existence of minimizers of 
our transport models. 


Proposition 2 . The discretized dynamic transport models and ( [Tol ) have a solution. 


Proof. We show that the proper, lower semicontinuous functions E and Ex are asympto¬ 
tically level stable which implies the existence of a minimizer. For the penalized problem, 
the asymptotic function Ex^oo reads 


^A,00 (^7 /) 


lim inf 

t—)-oo 


t 


We obtain 

t 


\ (||Jp(t(5^m',5f/' + i/+))||i + A||(D^|A) (^7') 


-||2 
2 


= \\Jp{Su.m',Sif + \n\\i + Xt\\{D^\Df) H 


tf Il2- 


Thus, (m, /) G ker(£’A,oo) implies 

(m, /) G ker(L)m| A), m G ker(5m), Sff > 0. 


( 12 ) 


For the constrained problem we have the same implications so that we can restrict 
our attention to the penalized one. By the definition of Sm we obtain ker(S'm) = 
[w ® 1 :w ^ for periodic boundary conditions and even N and ker(S'm) = { 0 } 


otherwise, where is defined as 1 = (1, —1,..., 1, —1)^ G In the case ker(*Sm) = {0}, 
the first and second condition in ( [T^ imply D^f = 0 so that by the definition of also 
/ = 0. In the other case, we obtain by the first condition in (12) that / = —DjD^m^ 
where = {DjDfY 




denotes the Moore-Penrose inverse of Df. Then 


Sff = —SfDlD^m = —SfDlD^{w 0 i) 

for some w G M^. Straightforward computation shows 

Sff = —SiD^D^ {w ®1 n) = —w (8) iw 

for some w G M^. Now the third condition in can only be fulfilled if w = 0 . 
Consequently we have in both cases 


*Sf/ = 0. 


(13) 
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Let p > 0, {/J^k}k be a bounded sequence and {{rrik^ fk)}k be a sequence fulfilling ©• 
By (12) and ([T^ we conclude 


Ex{{mk, fk) - p{rhj)) = \\Jp{SM,Sih + f+)\\i + X\\iD^m 

= Ex{{mk,fk))- 


Since {rrik^ fk) ^ lev(£’A5 fhis shows that (rrik^fk) — f) G as well 

and finishes the proof. □ 

Unfortunately, Jp{u^v) is not strictly convex on its domain as it can be deduced from 
the following proposition. 

Proposition 3. For any two minimizers (rui^fi), i = 1,2 of the relation 

Sih + f- ~ S,f2 + f- 


holds true. 

Proof. We use the perspective function notation from Remark For A G (0,1) and 
(ui^Vi) with Vi > 0, i = 1,2, we have (componentwise) 

Jp {X{ui,vi) + (1 - X){U 2 , V 2 )) = (Avi + (1 - A)?;2) ^ ) 

= (A^;i + (1 - A)^;2)V’ + a.iU-a% 2 f) 

and if ^ ^ ^ by the strict convexity of that 

Jp (A(«i,'yi) + (1 - X){U 2 ,V 2 )) < XJp{ui,Vi) + (1 - X)Jp{u 2 ,V 2 ). 

Setting Ui and Vi := Sffi + /“, i = 1, 2, we obtain the assertion. □ 

Remark 4. For periodie boundary eonditions, even N and fi = /o + yl, 7 G [0,min/o) 
the minimizer of @ is not unique. This ean he seen as follows: Obviously, we would 
have a minimizer (m, /) if m — w ®1 G ker(S'm) for some w G and there exists 
/ > 0 whieh fulfills the eonstraints. Setting := /(j — 1/2, A: = 0,..., P, these 

eonstraints read —2Pw ®1 = Thus, any w G sueh that 

= /o + 2wii, = /o + 2{wi + W 2 )l, p = fo + ‘2{wi + W 2 + ... + wp)i 

are nonnegative veetors provides a minimizer of We eonjeeture that the solution is 
unique in all other eases, but have no proof so far. 
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4. Primal-Dual Minimization Algorithm 

4.1. Algorithms 

For the minimization of our functionals we apply the primal-dual algorithm known as 
Chambolle-Pock algorithm mm in the form of Algorithm 8 in M- We use the 
following reformulation of the problems: 

Constrained Transport Problem: 

axgmm\\Jp{u,v)\\i +ic{m,f) (14) 

(nij) 

subject to SjnTTl = U, Sff + f~^ = V. 


Algorithm 1: Primal-Dual Algorithm for the Constrained Problem (14) 


Initialization: = 0^ = 0^ = bu^ = b^^ = 0, 0 G (0,1]^ 


r, a with ra < 1. 

Iteration: For r = 0,• • • iterate 


/^(r+l)\ 


argmin —1| 
imj)ec 2r 




-h ra 



2 . 

3. 

4. 


; ■“ 

^(r+l) 

^(r+1) 

r(r+l) 

r(r+l) 


argmin ||J,(u,^;)||i +^11 (“) 

+ y - 


f S'mTO('’+b\ _ f 0\ 

1 5f/(’-+b ) \fF 



Penalized Transport Problem: 

argmin||Jp(u,u)||i + A|| (DmlA) - /"111 (15) 

(m,/) '-v;- '\JJ 

A 

subject to S^m — S^f -\- = v. 

In the following we detail the first two steps of Algorithms and 

• Step 1 of Algorithm 1 requires the projection onto C, 

• Step 1 of Algorithm 2 results in the solution of a linear system of equations with 
coefficient matrix XA^A + y/ whose Schur complement can be computed via fast 
trigonometric transforms, 
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Algorithm 2: Primal-Dual Algorithm for the Penalized Problem © 

Initialization: = 0^ = 0^ = hu^ = = 0, 0 e (0,1]^ 

T, a with ra < 1. 

Iteration: For r = 0,• • • iterate 


/^(r+l)\ 


argmin^||(L>m|-Df) 

(mj) ^ 



/rii2 + ^i 




,siW\ 


2 — 4. as in Algorithmic 


• Step 2 of both algorithms is the proximal map of Jp. 


4.2. Projection onto C 


Step 1 of Algorithm 


1 requires to find the orthogonal projection of a := 

This means that we have to find a minimizer of \\Ax — f ~\\2 for 



which ||a — x ||2 attains its smallest value. Substituting y a — x we are looking for a 
minimizer y of ||— Aa +/“||2 with smallest norm \\y\\ 2 - By [HI Theorem 1.2.10], this 
minimizer is uniquely determined by (Aa — Therefore the projection of a onto C 
is given by 


nc(a) = a — A^ (^Aa — f ) (16) 

— a — A^(AA^)^ (Aa — /“) . 

Note that the projection onto C coincides with the one onto Cq if the given images /o 
and fi have the same mass. The Moore-Penrose inverse of the quadratic matrix AA^ is 
defined as follows: Let AA'^ have the spectral decomposition 

AA^ = Q diag(Aj) Q^. 


Then it holds 


(AA^)^ = Qdiag(Aj) with Xj := 


A if > if’ 

0 otherwise. 


The following proposition shows the form of (AA^)^^ in the one-dimensional spatial case. 
It appears that the projection onto C amounts to solve a two-dimensional Poisson equa¬ 
tion which can be realized depending on the boundary conditions by fast cosine and 
Fourier transforms in 0{NPlog{NP)) operations. 
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Proposition 5. Let Cn •= \/i (^j with cq := IjV^ and ej := 1, 

V \ / j,k=0 

j -j / —2nijk \ ^ 

j = 1,..., — 1 6e t/ie cosine matrix and Fn := W ( e n j be the N- 

th Fourier matrix. Set := and := (4sin^ • Then the Moore-Penrose 


inverse 


(AA^)^ in ([Ill) is given by 

{AA^)"^ I ^ ^^-i) diag(d) {Cp 0 Cn-i) mirror boundary, 

1 {Cp 0 Fn) diag(d) {Cp 0 F^) periodic boundary, 


where 


{ Ip ® A^^diag(djY^b\) + P^diag(dp^^’^) (g) In-i mirror boundary. 
Ip ® A^^diag(d^^) + P^diag(dp^’^^) ® In periodic boundary 

and dj := ^ if dj > 0 and dj — 0 otherwise. 

The proof is given in Appendix C. 


4.3. Schur Complement of XA^A + ^I 


To find the minimizer in Step 1 of Algorithm we set the gradient of the functional to 
zero which results in the solution of the linear system of equations 




XAy- + 



— ra 


(sifi:'’] 


Noting that XA^A+^I is a symmetric and positive definite matrix, this linear system can 
be solved using standard conjugate gradient methods. Alternatively, the next proposition 
shows how the inverse {XA^A + -I)~^ can be computed explicitly with the help of the 
Schur complement and fast sine,-, cosine- and Fourier transforms. The proposition refers 
to the one-dimensional spatial setting but can be generalized to the three-dimensional 
case in a straightforward way using the results of Appendix C. 

Proposition 6. Let Sn -1 — \f^ Then 

V V '' j:k=i ~ 

the inverse of the matrix XA^A F ^I is given by 

f I -X-^Y \ f X-i 0 \ / I 0\ 

\o I yVo I J ’ 

where 

i) for mirror boundary conditions 


Y = Df.®DN, 

X-^ =Ip® 5iv-idiag(AAr2d?^-L°i + p-'^AT-i, 

S~^ = {Sp-i (g) C'^)diag(^AP^dpi° (1 + rAA^^d™)“^ + ^fSp-i g Cn), 
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ii) for periodic boundary conditions 


Y = Dl® 

X-i = Ip® FNdmg{\N‘^d^j^^ + 

S-^ = (5p_i (8) Fjv)diag(Ap2dpi° ® (1 + rAiV^d^®'^)-^ + ® Fn). 

The proof is given in Appendix C. 


4.4. Proximal Map of Jp 

Step 2 of Algorithm ITl consists of an evaluation of the proximal map proxi of Jp. This 
can be done using the proximal map proxj* of J* and Moreau’s identity prox 0 (t) + 
prox 0 *(t) = t. Therefore, we state in the following first the dual function Jp. 

Lemma 7. For p G (1, +oc) and p + ^ = 1 have 


j;{a,b) 


0 if (a, b) G /Cp, 
+OC otherwise. 


where 

/Cp:=|(a,6)GM'^xM: J|a|'? + 6<o}. 

For the proof we refer to m or m Lemma 5.17]. After this preparation we are now 
able to compute proxij^. 

Proposition 8. Let p G (1, 2] and ^ ^ = 1. 

i) Then for x* G G M and a > 0 it holds 


proxi j 

a ^ 


(0,0) ii {ax*,ay*) eKp, 

otherwise, 


where 

h{z) := {ay* + 

and z G M>o is the unique solution of the equation 


z{l + h{z)) — cr|x*| = 0 (17) 

1 

in the interval [max( 0 , zq) ^,+oc), where zq := —qay"". 

ii) The Newton method eonverges for any starting point z > Zf) quadratieally to the largest 
zero of ( [TtI ) . 

Proof, i) By Moreau’s identity it holds 


prox0(t) + prox0*(t) = t 
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Figure 4: Projection onto the graph of the function (/)(x) = —^\x\^ for g = 3. 


and since we conclude 

{x,y) =proxi , (x*,^*) (18) 

(T ^ 

= (x - -prOX^j*(crX ,(77/ ), 

(7 P 

where proxij* = proxj* since is an indicator function. Note that by definition of 
Jp we have that ^ > 0 and ^ = 0 only if \x\ — 0. Now, proxj*((jx*,( jt/*) is the orthog¬ 
onal projection of ((jx*,(J7 /*) onto the set fCp (we could also compute the epigraphical 
projection of ((jx*,(J7 /*) onto the epigraph of (j){x) — ^\x\^ and reflect 7/, see also Fig- 
ureEj). If ((jx*,(j7/*) G /Cp, then proxj* ((jx*, (jt/*) = ((jx*,(J7 /*) and {x^y) — (0,0). So let 
((jx ,(J7/*) 0 JCp^ that means 


(77/* -h ^|(7X*|^ > 0. 


The tangent plane of the boundary of JCp in (x,7/) = (x,—^|x|^) is spanned by the 
vectors (e^, —|x|^“^x^)^, i = l,...,(i, where Ci G denotes the i-th canonical unit 
vector. Hence, the projection (x,7/) is determined by 7/ = —^|x|^ and 



= (7X* — X^ — ((77/* — 7/)|x|^ ^x^, i = l,...,(i 


so that 


(7X- 


Xi = 


-, i = l,...,d. 


1 + h{\x\)' 

Summing over the squares of the last equations gives 

^2 


I 12 I * 12 
X = X ' 


(7 


(i + ^(kl)) 


2 • 
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Since a solution has to fulfill 


y = y* -ly = y* + > o, 

it remains to search for the solutions with h{\x\) > 0. Now, h{z) > 0 is fulfilled for z > 0 
if and only if 


ay* + \z'i > 0, (19) 

1 

which is the case if and only if 2 ; > zq := max(0, —Then 2 ; := \x\ has to satisfy 
the equation 


2 ; (1 + h{z)) — (j|x* 


The function 


Lp{z) := z (1 + h{z)) — 


has exactly one zero in [ 2 : 9 ,+ 00 ). Indeed, by definition of 2:9 and (19) we see that 
V^(^o) < 0 , but on the other hand we have Lp{z) +00 as 2 ; ^ + 00 , so that ip has at 
least a zero in [ 2 : 9 , +oc). Since g > 2 for (p < 2 we have for 


1 


h'{z) = + {q- 2){ay* + -z^)z^-^ > 0 , 


2; > 2:9. 

Hence h and then also ip is strictly monotone increasing for 2 ; > 2 : 9 . Therefore ip has at 
most one zero in [ 2 : 9 ,+ 00 ). Finally, the assertion follows by plugging in {x^y) in (18). 

ii) Straightforward computation gives 


h”{z) = (Sq - 5)^2'^-^ + (q-2)iq- 3)(ay* + 

(p'{z) = 1 + h{z) + zh'{z), 

(p"{z) = 2h' {z) + zh"{z) 

= 3{q - l)z^^-^ + {q- l){q - 2){ay* + >0, ^ > zq- 

Since ip is monotone increasing and strictly convex for 2; > 2:9, the Newton method 
converges for any starting point 2; > 2:9 quadratically. □ 


5 . Numerical Results 


In the following we provide several numerical examples. In all cases we used P = 32 time 
steps and 2000 iterations in Algorithm and respectively. The parameters a and r 
were set to a = 50 and r = so that ar < 1, which guarantees the convergence of the 
algorithms. Of course, the algorithms do not use tensor products, but relations such as 
stated in (26) and their higher dimensional versions. The algorithms were implemented 
in Matlab and the computations were performed on a Dell computer with an Intel 
Core i7, 2.93 Ghz and 8 GB of RAM using Matlab 2014, Version 2014b on Ubuntu 
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Figure 5: Dynamic optimal transport between a red Gaussian and a blue one by the 
constrained model Q with different boundary conditions for the third (RGB) 
dimension. The initial and final images have the same mass. Top: mirror 
boundary conditions, bottom: periodic boundary conditions. 


14.04 LTS. Exemplary, the run time for 100 x 100 color images and 32 time steps varies 
between 10 and 15 minutes for the constrained and the penalized method, depending on 
the parameter choices for p, a, r and A. In our current implementation of the penalized 
method the fast transform approach is nearly as time consuming as the iterative solution 
of the linear system of equations with the CG method and an adequate initialization. 
If not explicitly stated otherwise, the results are displayed for p = 2, in which case the 
computation of the zeros in (0 slightly simplifies. 

With our first experiments we illustrate the difference between mirror and periodic 
boundary conditions in the color dimension, where at this point that the initial and 
the final images have the same mass. The images are displayed at intermediate time- 
points ^ where z = 0,..., 8. In Figure the transport of a red Gaussian into a 
blue one is shown, either with mirror or with periodic boundary conditions in the color 
dimension. Figure depicts the transport between two real images of polar lights. In 
order to have the equal mass constraint fulfilled, we first normalized both images to 
mass 1 and afterwards multiplied them with a common factor such that both images 
have realistic colors. Of course, this procedure works only if the initial and the final 
image share approximately the same mass. In both cases the use of periodic boundary 
conditions yields more realistic results. 

Further examples of the constrained model Q for several Gaussians and real images 
are given in Figures and i ^ The first row in Figures shows the transport of a red 
and a yellow Gaussian into a cyan and a blue Gaussian. The red and the blue Gaussian 
are spatially more extended compared to the yellow and the cyan one, but due to the 
fact that yellow and cyan have higher intensities, the masses of the red and blue Gaus¬ 
sians are approximately the same those of the yellow respective the cyan ones. This 
results in a very low interaction between the Gaussians during the transport, which is 
sightly visible in the background. Mainly, the red Gaussian is transported via a light 


^Images from Wikimedia Commons: AGOModra_aurora.jpg by Comenius University under CC BY SA 
3.0, Aurora-borealis_andoya.jpg by M. Buschmann under CC BY 3.0. 

^Images from Wikimedia Commons: Europe_satellite_orthographic.jpg and Earthlights_2002.jpg by 
NASA, Kohlbrandbrucke5478.jpg by G. Ries under CC BY SA 2.5, Kohlbrandbrucke.jpg by 
HafenCityl under CC BY 3.0. 
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Figure 6: Dynamic optimal transport between two polar lights by the constrained model 
Q with different boundary conditions for the third (RGB) dimension. The 
initial and final images have the same mass. Top: mirror boundary conditions, 
bottom: periodic boundary conditions. 


red into cyan and the yellow Gaussian is transported over violet to blue. The next row 
displays the intensity of the transported color images + G + B)^ in contrast to the 
(two-dimensional) transported intensity images in the third row. One sees slight differ¬ 
ences which arise due to the fact that the small mass difference can be transported only 
spatially and not through the color channels. 

The experiment is repeated in the fourth until sixth row, but this time the yellow and 
the cyan Gaussian are spatially more extended, thus having a significantly higher mass 
compared to the red respective the blue Gaussian. As a consequence, the interaction 
during the transport is higher, which is also clearly visible in the corresponding intensity 
images (fifth row). The color of the red Gaussian is transported similar as before, while 
the yellow Gaussian splits into two parts, one of them changing (as before) over violet 
to blue, while the other one goes over a light green to cyan. Further, as the Gaussians 
do not only travel in space, but also in color direction, the results are slightly smoother 
compared to the two-dimensional intensity transport, shown in the last row. 

Also for the real images, we assume that the images have the same mass. Indeed, the 
initial and final images had approximately the same overall sum of values, so that our 
normalization had no significant effect. In the first row of Figurea topographic map of 
Europe is transported into a satellite image of Europe at night. The second row displays 
the transport between two images of the Kohlbrandbriicke in Hamburg. In both cases 
one nicely sees a continuous change of color and shape during the transport. 

In Figure we give further examples for the transport of several Gaussians which may 
have different shapes in order to illustrate the transport of color and shape. Here, the 
initial and final images have different masses. The first row shows the transport of a 
yellow and a red Gaussian placed at the top of the image into a green and a blue Gaus¬ 
sian placed at the bottom. At this point, the yellow and the blue Gaussian are slightly 
more spatially extended compared to the red respective the green one. The red Gaussian 
changes over violet to blue. The yellow Gaussian, however, splits into two parts. While 
the main part is transported to green, a small part is separated and changes over orange 
to blue. In the second row, again a red and a yellow Gaussian are transported from the 
top into a green and a blue Gaussian at the bottom, but this time the yellow and the 
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Figure 7: Example for different color transitions obtained with the constrained model 
and periodic boundary conditions (first and fourth row), where the initial 
and final color images have the same mass. The second and fifth row show the 
corresponding intensity images, while the third and sixth row give the results 
obtained using two-dimensional transport of the initial and the final intensity 
images. 



Figure 8: Dynamic optimal transport between RGB images by the constrained model 
with periodic boundary conditions. The initial and final images have the 
same mass. 
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Figure 9: Example for different color and shape transitions by the constrained model © 
with periodic boundary conditions. The initial and final images do not have 
the same mass. 


green Gaussian are spatially more extended. Additionally, the Gaussians are no longer 
isotropic but have an ellipsoidal shape. In this case, additionally to the color also the 
shape changes continuously during time. Finally, the third row displays the transport 
of a white Gaussian into a yellow one. It shows that there appear no artificial colors 
during the transport, but the color transition proceeds as one may expect when looking 
at the RGB cube. 

Next we compare our approach with the approach of Rabin et al. m for microtex¬ 
tures, which is to the best of our knowledge the only approach that extends the dynamic 
optimal transport problem to a special class of color images. Note, however, that their 
approach is completely different from ours and works only for microtextures. At this 
point, microtextures are textures that fulfill the assumption of being robust towards 
phase randomization, in contrast to macrotextures, which usually contain periodic pat¬ 
terns with big visible elements (such as brick walls) or - more generally - the elements 
that form the texture-pattern are spatially arranged, see e.g. [29]. Based on the fact 
that microtextures can be well modeled as multivariate Gaussian distributions the au¬ 
thors of m propose to compute geodesics with respect to the Wasserstein distance 
W 2 between the Gaussian distributions that are estimated from the input textures /o 
and /i. This approach has the advantage that there exist closed-form solutions for the 
dynamic optimal transport between Gaussian measures. However, it is limited to the 
special class of microtextures, as natural images are not robust towards a randomization 
of their Fourier phase. In Figure [T0| we compare the results of our approach with the one 
for microtextures. In the case of microtextures both approaches yield similar results. 
Note that the approach of Rabin et al. m may fail for images which are orthogonal at 
some frequencies in Fourier domain. The second example demonstrates that the micro¬ 
texture technique m fails for natural images which possess contours and edges. 

Next, we turn to the penalized model (10). Figure 11 shows the influence of the reg¬ 
ularization parameter A when transporting a red Gaussian into a yellow one. Here, the 
initial and the final image have significantly different mass. The images are displayed at 










Figure 10: Comparison of our constrained model with periodic boundary conditions 
with the approach in EH for microtextures and the polar lights. In both 
cases, our RGB model is on top of the series and microtexture model is at 
the bottom. 


intermediate timepoints t i = 0,..., 8. The results change for increasing A from a 
nearly linear interpolation of the images to a transport of the mass. Further, for large A 
the results approach the one obtained with the constrained model Q, which is reason¬ 
able. 

Finally, we consider the influence of the parameter p G (1,2]. The corresponding results 


for our penalized model (10) are given in Figures |12| and |13[ In all our experiments we 
observed only rather small differences. Note however that in m Wasserstein barycen- 
ters were considered for p = 1, 2, 3 which show significant differences. 

Further examples and videos, in particular for real images, can be found on our website 
http://www.mathematik.uni-kl.de/imagepro/members/laus/color-OT, 


6. Summary and Conclusions 

Our contribution can be summarized as follows: 

i) We propose two discrete variational models for the interpolation of RGB color 
images based on the dynamic optimal transport approach. To this end, we consider 
color images as three-dimensional objects, where the “RGB direction” is handled 
in a periodic way. We focus on a discrete matrix-vector approach. 

ii) Our first model relaxes the continuity constraint so that a transport between images 
of different mass is possible, while the second model allows even more flexibility by 
just penalizing the continuity constraint with different regularization parameters. 
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Figure 11: Comparison of penalized and constrained color optimal transport (from top to 
bottom): penalized optimal transport for different regularization parameters 
A G {0.1,1,10,100} and constrained optimal transport. 



Figure 12: Dynamic optimal transport of RGB images using the penalized model (10) 
with periodic boundary conditions. Comparison of p = 2 (top) and p = 1.5 
(bottom) for A = 1. The images are displayed at intermediate timepoints 
t= |,i = 0,...,4. 


24 













Figure 13: Dynamic optimal transport of RGB images using the penalized model (10) 
with periodic boundary conditions and A = 1 for p = 2 and p = 1.5. From 
left to right: Initial images /o, /i, result for p = 2 and p = 1.5 at time t = 0.5 
and the absolute difference between the two results. 


hi) We provided an existence proof and a brief discussion on the uniqueness of the 
minimizer. 

iv) Interestingly, the step in the chosen primal-dual algorithm which takes the continu¬ 
ity constraint into account requires the solution of four-dimensional Poisson equa¬ 
tions with simultaneous mirror/periodic (constrained model) or zero/mirror/periodic 
boundary conditions (penalized model). Here, fast sine, cosine and Fourier trans¬ 
forms come into the play. 

v) We consider the case p G (1,2] and give a careful analysis of the proximal mapping 
of Jp, p G (1,2]. This includes the determination of a starting point for the 
Newton algorithm to ensure its quadratic convergence and a stable performance 
of the overall algorithm. 

vi) We show numerous numerical examples. 

There are several directions for future work. 


1) One possibility is to add several additional priors. So far, the present model has 
difficulties to transfer sharp contours. A remedy could be the penalization of a 
total variation (TV) term with respect to /, which results for some 7 > 0 in the 
functional 

argmin{||Jp(S'mm,S'f/ + /+)||i + 7 TV(/)}. (20) 

{mj)ec 

In the following we use a spatial TV term, summed over time, i.e., in one dimension 

TV(/) = ||(/p®Z)^)/||i. 


Figure shows the performance of such a model in the one-dimensional case. 
Here, the approach without the TV term leads to some overshooting and blurring 
at the edges. With the TV term the transport takes place in a more natural way, 
in particular the sharp edges are preserved. For the one-dimensional example this 
works well. However, in higher dimensions one has to be more careful. In Fig¬ 


ure 15 a comparison of an isotropic and an anisotropic TV regularizer is shown. 
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Figure 14: Transport of a one-dimensional signal with sharp edges using a TV penalized 
functional ( [ 20 | ) (left), where 7 = 0.03 and the result without TV regularization 
(right). 


The isotropic one leads, as could be expected, to a rounding of the corners. The 
anisotropic regularizer prefers horizontal and vertical edges. In this way, the shape 
of the object is preserved during the transport. Note that due to the smeared 
boundary the square appears to be smaller. To preserve the shape of arbitrary 
transported objects, one would have to adjust the regularizer according to the di¬ 
rection of the edges. 

The idea of penalizing TV terms for the transport can be found for gray-value im¬ 
ages, e.g. in [I21I331. For image denoising a Wasserstein-TV model was successfully 
applied in [131 El M • 

2) Using a barycentric approach the interpolation of microtextures in [41] works also 
between more than two images. So far this task can not be handled via the 
dynamic optimal transport approach. One idea might it be to formulate a dynamic 
barycenter optimal transport problem or to use a multimarginal model, see e.g. 
Section 1.4.7. in |43| and the references therein. 

Acknowledgement: Funding by the DFG within the Research Training Group 1932 
is gratefully acknowledged. 


A. Diagonalization of Structured Matrices 


In the following we collect known facts on the eigenvalue decomposition of various differ¬ 
ence matrices. For further information we refer, e.g., to [39l[l6|. The following matrices 
Cn and Sn are unitary, resp., orthogonal matrices. The Fourier matrix 


Fn- = 



— 27Tijk 

e ^ 


n 

j,k=0 
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Figure 15: Transport of a two-dimensional square with sharp edges using a TV penalized 
functional, where 7 = 0.05. First row: isotropic TV, second row: anisotropic 
TV, third row: no TV regularization. 


diagonalizes circulant matrices, i.e., for a := G we have 


/ «0 

dn—l 

. . di 

\ 


ai 

do 

.. a2 







= Fndmg{^/nFna)Fn 

\ ^n—1 

di 

.. do 

) 



— -^ndiag( 


In particular it holds 

Ar := ^{D^rnr = ^oriorr 


( ' 

-1 2 -1 

1-1 

1 

1-1 

1 

-1 2-1 

-1 2 / 


= ^ndiag(dP®'')F„ 


( 21 ) 


( 22 ) 


with := (4sin^^)7o- The operator typically appears when solving the 

one-dimensional Poisson equation with periodic boundary conditions by finite difference 
methods. 







The DST-I matrix 


and the DCT-II matrix 


:= 


. jkiT 


ej cos 


j{2k + l)7r 


with eo := and Cj := 1, j = 1 ,..., n — 1 are related by 

Dn = -Sn-I (o I diag(d^®!i)^^ Cn, (23) 

where := (4 sin^ Further they diagonalize sums of certain symmetric 

Toeplitz and persymmetric Hankel matrices. In particular it holds 


A zero _7^ t^t 
^n-1 — 


2 -1 

-1 2 -1 


-1 2 -1 
-1 2 


= 5^_idiag(d^^4i)5n-i 


A mirr ijT ij 


1 -1 

-1 2 -1 


-1 2 -1 
-1 1 


= C^diag(dr^)a 


with d™^^^ := ( ) = (dsin^ ) . The operators and are related to 

V^n-iy ^ 

the Poisson equation with zero boundary conditions and mirror boundary conditions, 
respectively. 


B. Computation with Tensor Products 

The tensor product (Kronecker product) of matrices 


and B = 


28 



is defined by 


A®B 


ai,iB 


^l,nB \ 




The tensor product is associative and distributive with respect to the addition of matri¬ 
ces. 

Lemma 9 (Properties of Tensor Products). 

i) {A (g) By = for A G B G 

Let A^C ^ £m,m B^D ^ Then the following holds: 

ii) {A ® B){C ® D) = AC ® BD for A,C e and B,D e C^’^. 

iii) If A and B are invertible, then A® B is also invertible and 

{A^B)-^ = A-^ ^B-^. 


The tensor product is needed to establish the connection between images and their 
vectorized versions, i.e., we consider images F G columnwise reshaped as 

/:=vec(F) G 

Then the following relation holds true: 

Yec{AFBy = {B® A)f. (26) 


C. Proofs and Generalization of the Tensor Product Approach 
to 3D 


Proof of Proposition By definition of A and using (25), (22), we obtain for periodic 
boundary conditions 


AA^ = Ip^ + Dj^Dp ® In 

= Ip0 (g, /jv 

= {CJ> ® Fn) diag(/p ® ^ ^ 


Similarly we get with (25) for mirror boundary conditions 


AA^ = Ip(g> DI_^Dn-i + D},Dp ® In-1 
= Ip(S) ® In-1 

= (CJ, ® C^_i) diag(Ip ® N^dtL\ + P‘^d^L\ ® In-i) {Cp ® Civ-i) 
which finishes the proof. 


□ 
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Proof of Proposition^^ By definition of A we obtain 

XD^D^ + n XD^Di 


+ -I = 


T V XDJDf + ^I 

SO that the inverse can be written by the help of the Schur complement 

S :=Z-Y'^X-^Y 


T n. I 1 


X Y 

yx y 


as 


X Y 


-1 


I -X-^Y 
0 I 


x-i 

0 


0 

5-1 


I 

_yTx-i 


By (22) and ( [2^ we have with D G that 

X-i = {XDlD^ + = {Ip ® XN^DD^ + p)"^ 

= /p ® {XN‘^DD^ + p)-i 

Ip ® S'w-idiag(AA^^d^f°]^ + ^)~^Sn-i mirror boundary, 

Ip (g) F/vdiag(AA^^d^^ + periodic boundary. 

The Schur complement reads as 

5 = {XDjDf + p) - X^DJD^X-^DID{ 

= {XDpDl ®In + p) - X^{Dp ® D^) {Ip 0 {XN^DD'^ + p)-i) {D}> 0 D) 
= {XDpDj, (g) liv + p) - X^{DpD}> (g) D'^{XDD'^ + piv)"^y>) 

= XDpDj, ® {In - XD'^{XDD'^ + ^In)~^D) + p. 

By ( [2T| ) we have 


{D^^^y = NFNdmg{-l + e+^^^’^/^)kFN 


and 


DP"’’ = XFjvdiag(-l + 

SO that we obtain for periodic boundary boundary conditions 

Ik - MD^r(\D’^’{D%‘-r + \In)-^D'S‘ = F„diag (1 + rj^d^Y ^N- 

Therewith it follows with ( [24l ) 

S = 5p_idiag(AF^dpi° )5p_i (g) Fjvdiag (l + Fn + p 

= (5p_i ® FAr)diag(Ap 2 dpi°i ® (1 + + i)(5p_i ® P^v) 

which yields the assertion for S~^ in the periodic case. 

For mirror boundary conditions we compute using (23) 

5 = (5p_i (g C^)diag(Ap2dF-i ® (1 + rP^dr'')-! + P(5p_i (g Cn) 

and inverting this matrix finishes the proof. 


□ 
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Discretization for three spatial dimensions + time: For RGB images of size 
X A^2 X ^ 3 , where — 3, we have to work in three spatial dimensions. Setting 
N := (A^i, ^ 3)5 j {ji^j2^j3) and defining the quotient componentwise we obtain 


fi = (/*e^)). ,,,, e ^ = 0 , 1 , 




/ = (/("^ 


A 
’ p 


N,P-1 




• m = (mi, m 2 , m 3 ), with 



J2-1/2 

N 2 


J3-1/2 

3 


k-1/2 

P 


Wi-1,W2,3,P 
J1=1J2 = 1,J3 = 1A=1 


G M^i-bA^ 2 , 3 A^ 


Ji- 1/2 J 3 - 1/2 k-1/2 

Wi ’ W2 ’ 3 ’ P 


Wi,W2-1,3,P 
J1 =1J2 = 1,J3 = 1A=1 




Ji- 1/2 

Ni 


J2-1/2 

W 2 


J3 fe-1/2 
3 ’ P 


Wi,W 2 , 2 ,P 

jl = lj 2 = lj3=0A=l 


^ j^Wi,W2,3,P^ 


In the definition of m we take the periodic boundary for the third spatial direction into 
account. Analogously as in the one-dimensional case, when reshaping m and / into 
long vectors, the interpolation and differentiation operators can be written using tensor 
products. For the interpolation operator we have 

/ {Ip (g) /3 (g) In2 ® 

Smm = (/p (g) /3 ® S'w 2 ® ^iVi)P ^2 

\(/p (g) 53 (g) In2 ® /Ari)p^3/ 


and 


Sff = {Sp 0 / 3 ® In2 ® lNi)f, 

which means, that computes the average of m^ with respect to the i-th coordinate, 
i = 1, 2, 3, and Sff computes the average of / with respect to the time variable. Similarly 
we generalize the difference operator. Then, reordering / and m into large vectors, the 
matrix form of the operator A is 

A = {ip 0 / 3 ® In2 ® I /p ® ^3 ® ® /wi I 

Ip ® (Df ® In2 ® I Dp ® /3 ® 7^2 ® Ini) 


so that reads as 

AA^ = /p ® 73 ® 7^2 ® + 7p ® 73 ® ® 7^1 

+ 7p ® (7^3 (7^3 ® 7^2 ® Ini + T^pT^p ® 73 ® 7^2 ® 7wi 

= ® F3 ® C^^_i ® C^^_i) diag(d) (Cp ® F3 ® CW2-1 ® Cn,-i) , 
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where 


d l3P(N2-l) ® iVfdiag(d™1i) + hp (g) iV2^diag(d™i) (g) Ini-i 

+ Ip® 3^diag(d|®") (g) /(7V2-i)(jVi-i) + P^diag(d^'"") (g) l3P(7V2-i)(jVi-i). 

For the three-dimensional spatial setting we have to solve a four-dimensional Poisson 
equation, which can be handled separately in each dimension. For the constrained prob¬ 
lem, this can be computed using fast cosine and Fourier transforms with a complexity 
of 0 {NiN 2 Plog{NiN 2 P)). 
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